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Synchronization of coupled oscillators is fre- 
quently observed in nature and technology 1 ' 2 . 
Recently, the study of synchronization phenom- 
ena in complex networks has received much 
attention 3 ) 4 !^! 7 ^?! 1 '?! 11 ! 12 . A classical model for 
the phase dynamics of weakly coupled oscillators 
is that of Kuramotc 13-14 , who showed that as the 
coupling strength is increased there is a transi- 
tion from incoherent behavior to synchronization. 
The Kuramoto model assumes all-to-all connec- 
tivity and positive coupling (i.e., the coupling of 
two oscillators tends to reduce their phase dif- 
ference). However, it has been recently noted 
that the topology of real world networks is of- 
ten very complex. In the current paper, gener- 
alizing our previous work which considered the 
case of undirected coupling networks with posi- 
tive coupling 12 , we discuss the synchronization of 
phase oscillators interacting on directed networks 
with mixed positive/negative coupling. 

I. INTRODUCTION 

The classical Kuramoto modeli&ii describes a collec- 
tion of globally coupled phase oscillators that exhibits 
a transition from incoherence to synchronization as the 
coupling strength is increased past a critical value. Since 
real world networks typically have a more complex struc- 
ture than all-to-all couplingi^*i&, it is natural to ask what 
effect interaction structure has on the synchronization 
transition. In Ref.0, we studied the Kuramoto model al- 
lowing general connectivity of the nodes, and found that 
for a large class of networks there is still a transition to 
global synchrony as the coupling strength exceeds a criti- 
cal value k c . We found that the critical coupling strength 
depends on the largest eigenvalue of the adjacency ma- 
trix A describing the network connectivity. We also de- 
veloped several approximations describing the behavior 
of an order parameter measuring the coherence past the 
transition. This past work was restricted to the case in 



which A nm = A mn > 0, that is, undirected networks in 
which the coupling tends to reduce the phase difference 
of the oscillators. 

Most networks considered in applications are 
directedi^i^, which implies an asymmetric adja- 
cency matrix, A nm ^ A mn . Also, in some cases the 
coupling between two oscillators might drive them to be 
out of phase, which can be represented by allowing the 
coupling term between these oscillators to be negative, 
A nm < 0. The effect that the presence of directed 
and mixed positive/negative connections can have on 
synchronization is, therefore, of interest. Here we show 
how our previous theory can be generalized to account 
for these two factors. We study examples in which cither 
the asymmetry of the adjacency matrix or the effect 
of the negative connections are particularly severe and 
compare our theoretical approximations with numerical 
solutions. 



II. BACKGROUND 

In Ref. ^3 we considered the onset of synchronization 
in networks of heterogeneous coupled phase oscillators. 
This situation can be modeled by the equation, 

N 

On =UJ n + k'^2 A nm Sin(9 m ~ 6 n ), (1) 
m— 1 

where 9 n , uj n are the phase and natural frequency of oscil- 
lator n, and N ^> 1 is the total number of oscillators. The 
frequencies uj n are assumed to be independently drawn 
from a probability distribution characterized by a den- 
sity function g(uj) that is symmetric about a single lo- 
cal maximum at to = ZJ. The mean frequency ZJ can be 
shifted to ZJ = by introduction of the change of variables 
On — ¥ n —ZJt. Thus we henceforth take ZJ = 0. The adja- 
cency matrix {A nm } determines the network connecting 
the oscillators. Positive coupling was imposed in Ref. 12 
by the condition A nm > 0. Furthermore, the matrix A 
was assumed to be symmetric and thus only undirected 
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networks were considered. In this Section we will review 
our results for this class of networks, following Sec. II of 
Ref. IT^ Thus throughout this Section A nm — A mn > 0. 

In order to quantify the coherence of the inputs to a 
given node, a positive real valued local order parameter 
r n is defined by 



N 



(2) 



where (...)* denotes a time average. To characterize the 
macroscopic coherence for the whole network, a global 
order parameter is defined by 



T N r 

l^n=l 



where d n is the degree of node n defined by 

N 



(3) 



dn ^ A nm . 



m— 1 



In terms of r„, Eq. Q can be rewritten as 
9 n = u n - kr n sin(#„ - tf) n ) - kh n (t), 



(4) 



(5) 



where the term h n (t) takes into account 
time fluctuations and is given by h n = 
Im{e- ie "J2m A nm{(e l9m }t-e ie ™)}, where 7m stands 
for the imaginary part. As argued in Ref. IT^ , when the 
number of connections into each node is large, the term 
h n is small compared to r n and we obtain approximately 



(see Appendix A of Ref. IT^ . and we obtain from the real 
and imaginary parts of Eq. (JSJ 



A nm cos(ip m - ip n )\ l 



\<kr„ 



kr r . 



(9) 



^ A mn sin(V> m - tpn) ' 



and 



\uj m <kr 
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\u) m \<krr, 



A nm COs(lp m -i>n) ( ) (10) 



^2 A nm sin(^„ 

\uJ m \<kr m 




^>m 
kTr. 



Introducing the assumption that the solutions ip n , "fn are 
statistically independent of u n (see Ref. ll2|) and using the 
assumed symmetry of the frequency distribution g(u>) we 
obtain from Eq. JSJ the approximation, 



^ ^ Anm 

, n \<kr m 



COs(lpm - Ipn 



kTn 



, (11) 



and the right side of Eq. (|1U[) is approximately zero for 
large number of connections per node. The solution of 
Eq. (|llfl with %p n = ip m for all n is the one corresponding 
to the smallest value of k, and thus corresponds to the 
smallest critical coupling k c leading to a transition to 
a macroscopic value of r n . Therefore we consider the 
equation 



w„ - kr n sin(6»„ - if> n )- 



(6) 



Henceforth, we will assume that the number of connec- 
tions into each node is large enough that we can neglect 
the time fluctuations represented by the term h n - For a 
discussion of the effect of nodes with few connections, see 
Sec. VI of Ref. [H 

From Eq. ©, we conclude that oscillators with \cu n \ < 
kr n become locked, i.e., for these oscillators 9 n settles at 
a value for which 



sin(6> n - ip n ) = uj n /(kr n ). 



Then 



E ^ 



-V>») 



(7) 



(8) 



I </cr T) 



\ujjn \ >kr m 

The sum over the non-locked oscillators can be shown to 
vanish in the large number of connections per node limit 



m.\<kr„ 



A, 



We refer to this approximation [Eq. (|12|l ]. based on ne- 
glecting the time fluctuations, as the time averaged theory 
(TAT). In Ref. [l2| we showed numerically that this ap- 
proximation consistently describes the large time behav- 
ior of the order parameter r past the transition for vari- 
ous undirected networks with positive coupling strengths 
(i.e., A nm = A mn > 0). 

Averaging over the frequencies, one obtains the fre- 
quency distribution approximation (FDA): 



r„ = k 



g(zkr m )\/l ~ z 2 dz. 



(13) 



The value of the critical coupling strength can be ob- 
tained from the frequency distribution approximation by 
letting r n — > + , producing 



„(o) _ ±\^ A 



, m 



(14) 
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where fco = 2/[7rg(0)]. The critical coupling strength thus 
corresponds to 



k 

A ' 



(15) 



where A is the largest eigenvalue of the adjacency matrix 
A and r^ is proportional to the corresponding eigenvec- 
tor of A. By considering perturbations from the crit- 
ical values as r„ = r n ^ + 8r n , expanding g(zkr m ) in 
Eq. (|13fl to second order for small argument, multiply- 
ing Eq. |Q3J by r { n } and summing over n, we obtained 
an expression for the order parameter past the transition 
valid for networks with relatively homogeneous degree 
distributions^: 



Vi 



- 1 



for < (k/k c ) — 1 <C 1, where 



(u) 2 X 2 



N(d) 2 (u 4 



(16) 



(17) 



a = — ng"(0)ko/16, u is the normalized eigenvector of 
A corresponding to A, and (...) is defined by (x q ) — 

The mean field theory (MFT)2iiS was obtained from 
the frequency distribution equation by introducing the 
extra assumption that the local mean field is approxi- 
mately proportional to the degree, r n — rd n . Substitut- 
ing this into Eq. I|13[l and summing over n we obtained 



iV iV 1 

Y / d m = kY / d lrl g(zkrd m )^l-z 2 dz. (18) 

m—1 m— 1 1 

Letting r — » + , the critical coupling strength is given by 



An expansion to second order yields 
% \ ( k 



r 2 = 



- 1 



v / \ A ™/ 

for < (k/k m f) — 1 <C 1, where 

(d 2 ) 3 



k m f 



1)2 



(19) 



(20) 



(21) 



Comparing the above three approximations, we note 
the following points: 

1. The TAT requires knowledge of the adjacency ma- 
trix and the particular realization of the oscillator 
frequencies w n at each node. 



2. The FDA requires knowledge of the adjacency ma- 
trix and the frequency distribution, but averages 
over realizations of the node frequencies. 

3. The MFT (like the FDA) averages over realizations 
of the node frequencies, but only requires knowl- 
edge of the degree distribution d m (knowledge of 
the adjacency matrix is not required). 

4. Computationally, the TAT and the FDA are more 
demanding than the MFT; all three, however, are 
much less costly than direct integration of Eq. ||TJ 
to find the time asymptotic result. 

5. Finally, one might suspect that the TAT is more ac- 
curate for describing a specific system realization, 
given that one has knowledge of the network and 
the realization of the oscillator frequencies oj n on 
each node, while the FDA might be more appro- 
priate for investigating the mean behavior averaged 
over an ensemble of realizations of the oscillator fre- 
quencies. 



III. DIRECTED NETWORKS 

In this Section we will extend our previous results to in- 
clude directed networks, A nm ^ A mn . As in the previous 
Section, we will assume that the number of connections 
per node (both incoming and outgoing) is large, that the 
frequencies are drawn randomly from a distribution sym- 
metric around its unique local maximum at lo — 0, and 
that the coupling is positive, A nm > 0. We define the 
in- degree d™ and out- degree of node n as 



and 



N 

^ ^ Anm 

m—1 
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(22) 



(23) 



For directed networks, the degrees and d° n ut may be 
unequal, and it is therefore necessary to take this dif- 
ference into account when developing approximations for 
the synchronization transition based on the degree of the 
nodes [e.g., the mean field theory, Eq. lfT5)l]. 

The approximations to r given by the time averaged 
theory [Eq. lfP2)l]. the frequency distribution approxima- 
tion [Eq. (|13fl ], and the estimate for the critical coupling 
constant given by Eq. I|15f) ar e still valid in this more 
general case. The existence of a nonnegative real eigen- 
value A larger than the magnitude of any other eigen- 
value is guaranteed for matrices with nonnegative entries 
by the Frobenius theorem^, and we use this eigenvalue 
in Eq. (Q2J). 

We now consider the perturbation solution to the FDA 
[Eq. I|13|) ] for (k — k c ) small taking into account asym- 
metry of A. Expanding Eq. (|13|) to second order in kr n , 
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inserting r n — r n ^ + 5r n , and canceling terms of order 
, the leading order terms remaining are 



Sr„ 



' Y. A nm5r m -^rY. A ^f (24) 



+ - 



1 



k r X 



r (0) 
nm ' m ' 



In order for Eq. I|24l) to have a solution for 8r n , it must 
satisfy a solubility condition. This condition can be ob- 
tained as follows. Let u n be an eigenvector of the trans- 
pose of A, A T , with eigenvalue A. Multiplying Eq. I|24|) 
by u n , summing over n and using Eq. I|14|l , we obtain 



X/ m ( r ™ )^u m k k c 



(0) 



\k :i 



(25) 



In terms of u and u, eigenvectors of A and A T associated 
with the eigenvalue A, the square of the order parameter 
r can be expressed as [cf. Eqs. (|To)l and l(T7l) ] 



Vi 



-3 



for < (k/k c ) — 1 <C 1, where 

(u) 2 (uu)X 2 



N(d)< 



(26) 



(27) 



and (x p y q ) is defined by (x p y q ) = J2n=i x nVn/ N - We wil1 
refer to this generalization of the perturbation theory as 
the directed perturbation theory (DPT). 

The mean field theory can also be generalized for di- 
rected networks by introducing the assumption r n = 
rd™ . We obtain as a generalization of Eq. H18fl the di- 
rected mean field theory (DMFT) 



N N l 

£C = *E /_ 9(zkratWl-z 2 dz. (28) 

m—1 rn— 1 — 

Letting r — > + , the critical coupling strength is given by 



k — k m f — ko 



(d in d° 



(29) 



An expansion to second order yields [cf. Equations 120|) 
and 



{ ak 2 ) \k m f J \k m f 



for < (k/k m f) — 1 <C 1, where 



in Jout\3 
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((d m ) 3 d° ut }(d in ) 2 ' 



(30) 



(31) 



IV. NETWORKS WITH NEGATIVE COUPLING 

Here we extend our previous results to the case in 
which the matrix elements A nm are allowed to be nega- 
tive. In this solution to Eqs. I© and HlOfl in which 
all the phases are equal, (ip n = ip m for all n,m), does not 
necessarily exist. [In fact, if one were to set ip n = ip m in 
Eq. the right hand side of Eq. I|12|l could be negative, 
while by definition r n is nonnegative.] 

Although in this section we will assume k > 0, the 
case k < can be treated by redefining k — > — k and 
A nm — > —A nm - By neglecting the contribution of the 
drifting oscillators, using the symmetry of g(uj) and the 
assumed independence of ip n and r n from uj n , we obtain 
from Eqs. (0), J7J) and © the equation 



= E ^ 



<fcr 




(32) 



Our approach will now be to solve Eq. 1)32(1 numerically 
for ip n and r„ . We note that such numerical solution will 
still be orders of magnitude faster than finding the exact 
temporal evolution of the network by numerically inte- 
grating Eqs. (Q. In order to numerically solve Eq. i|32|) 
for the variables tp ni r„, we look for fixed points of the 



following mapping, (r n ,tp n ) 



4> n +1 ), defined by 



rl+'e 1 ^ 1 




I in 



Repeatedly iterating the above map starting from ran- 
dom initial conditions, the desired solution will be pro- 
duced if the orbit converges to a fixed point. We will 
discuss the convergence of this procedure when consider- 
ing particular examples. 

We now comment on some aspects introduced by con- 
nections with negative coupling. First, we note that when 
the coupling between the oscillators is positive, the ef- 
fect of the coupling between them is a tendency to re- 
duce their phase difference. In this case, as k — > oo, the 
phases synchronize, 9 n — > 0. There is in this case fre- 
quency and phase synchronization [i.e., -^(dn — m ) — > 
and (0 n — 8 m ) — * 0]. On the other hand, two oscillators 
coupled with a negative connection A nm < tend to os- 
cillate out of phase. However, in a network with many 
nodes and mixed positive/negative connections, the rel- 
ative phases of two oscillators can not in general be de- 
termined only from the sign of their coupling. When the 
oscillators lock, their relative phase is determined by ipn 
[let k — ► oo in Eq. Q], and in general the phases ip n can 
be broadly distributed in [0, 2n). Therefore in this case 
we expect frequency synchronization, but not phase syn- 
chronization [i.e., -§:{0 n — 9 m ) — > but (9 n — 9 m ) 0]. 
We also note that in this case the order parameter r, as 
we have defined it in Eq. @, may attain values higher 
than 1 for k — > oo. We therefore replace the definition 
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© by 



T N r 



(34) 



Note that if ^4„ m > this definition reduces to the pre- 
vious one. 

From Eq. we have for k — > oo 



E m ,„ ^nm COs(lp m - Ipn) 
Em,n l-^nm I 



(35) 



The order parameter achieves its maximum value, r = 1, 
when the phase difference ip m — ip n between two oscilla- 
tors is for positive coupling (A nm > 0) and ir for neg- 
ative coupling (A nm < 0). An order parameter smaller 
than 1 as k — > oo indicates frustration in the collection 
of coupled oscillators, i.e., the phase difference favored 
by the coupling between each pair of oscillators cannot 
be satisfied simultaneously by all pairsiS. The order pa- 
rameter is similar to the overlap function used in neural 
networks for measuring the closeness of the state of the 
network to a memorized pattern^. 

Using the assumption that the number of connections 
per node is large, we average Eq. over the frequencies 
to obtain the approximation 

N i 

r n e**» = k V A nm e l ^r m / y/l - z 2 g{zkr m )dz.{36) 

m=l J- 1 

The critical coupling strength k c can be estimated by 
letting r n — > + to be as in Sec. [H] 



ko 
A ' 



(37) 



where fco — 2/[7rg(0)] and we have assumed the existence 
of a positive real eigenvalue A which is larger than the 
real part of all other (possibly complex) eigenvalues of 
A. We now discuss the validity of this assumption. 

If the adjacency matrix A is asymmetric and there are 
mixed positive/negative connections (both A nm > and 
A n ' m > < for some n,m,n' ,m'), it might occur that the 
matrix A has no real eigenvalues, or it has complex eigen- 
values with real part larger than the largest real eigen- 
value. In our examples we find, however, that when there 
is a bias towards positive coupling strengths, there is a 
real eigenvalue A with real part larger than that of the 
other eigenvalues. Furthermore, the largest real part of 
the remaining eigenvalues is typically well separated from 
A. This issue is discussed further and illustrated with the 
spectrum of a particular matrix in Appendix lAl 

So far, we have considered situations in which coupling 
from oscillator m to oscillator n favors a phase difference 
Qn — 6m = (positive coupling, A nm > 0), or situations 
in which a phase difference 9 n — 8 m = ir is favored (neg- 
ative coupling, A nm < 0). A more general case is that 
in which coupling from oscillator m to oscillator n favors 
a phase difference n m — oz nm , with ^ oz nm < 2tt. 



(Such nontrivial phase differences could be favored, for 
example, by a time delay in the interaction of the oscil- 
lators in conditions in which, in the absence of a delay, 
their interaction would reduce their phase difference to 
zero.) This more general case can be described by the 
following generalization of Eq. : 



N 



k \ A r. 



+ a nm ). (38) 



In this scenario, positive coupling corresponds to a nm — 
and negative coupling to a nm = n. By considering 
complex values of the coupling constants, 



Anm — | A nrn | e 



(39) 



the same process described at the beginning of this Sec- 
tion can be used to show that Eq. I)32fl is still valid in 
this more general case. For simplicity, in our examples 
we will consider cases in which a nm is either or tt. 



V. EXAMPLES 

In this section we will numerically test our approxima- 
tions (Sees. 11111 and HV|) with examples. 

In Ref. [13 we showed how our theory described the 
behavior of the order parameter r for a particular re- 
alization of the network and the frequencies. Although 
the agreement was very good, there was a small but no- 
ticeable difference between the time averaged theory and 
the frequency distribution approximation. Here, besides 
the asymmetry of the adjacency matrix, we will investi- 
gate the variations that occur when different realizations 
of the network and the frequencies of the individual os- 
cillators are considered. We will show that the small 
discrepancies mentioned above can be accounted for by 
averaging over many realizations of the frequencies. 

We will compare the approximations described in this 
section with the numerical solution of Eq. for differ- 
ent types of networks. When numerically solving Eq. 
the initial conditions for 9 n are chosen randomly in the 
interval [0, 2-7r) and Eq. is integrated forward in time 
until a stationary state is reached (stationary state here 
means stationary in a statistical sense; i.e., although the 
solution might be time dependent, its statistical proper- 
ties remain constant in time). From the values of 9 n (t) 
obtained for a given k, the order parameter r is estimated 
using Eqs. @ and ©, where the time average is taken 
after the system reaches the stationary state. (Close to 
the transition, the time needed to reach the stationary 
state is very long, so that it is difficult to estimate the 
real value of r. This problem also exists in the classi- 
cal Kuramoto all-to-all model.) The value of k is then 
increased and the system is allowed to relax to a sta- 
tionary state, and the process is repeated for increasing 
values of k. Throughout this section, the frequency dis- 
tribution is taken to be g(u>) = |(1 — lo 2 ) for \uj\ < 1 and 
otherwise. 



6 



A. Example (i), A Randomly Asymmetric Network 
with A > 

As our first example [example (i)] we consider a di- 
rected random network generated as follows. Starting 
with N ^> 1 nodes, we consider all possible ordered pairs 
of nodes (n, m) with n ^ m and add a directed link 
from node n to node m with probability s. (Equivalently, 
each nondiagonal entry of the adjacency matrix is inde- 
pendently chosen to be 1 with probability s and with 
probability 1 — s, and the diagonal elements are set to 
zero.) Even though the network constructed in this way 
is directed, for most nodes « e?° u *. For N = 1500 

„2 1-0 




10 realization!; 



[ realization 



k/k f 



1.6 



FIG. 1: (a) Average of the order parameter r 2 obtained from 
numerical solution of Eq. Q over 10 realizations of the net- 
work and frequencies (triangles) , from the frequency distribu- 
tion approximation (solid line) and from the directed mean 
field theory (long dashed line) as a function of k/k c . (b) Or- 
der parameter r 2 obtained from numerical solution of Eq. Q 
for a particular realization of the network and frequencies 
(boxes), from the time averaged theory (short dashed line) 
and from the frequency distribution approximation (solid line) 
as a function of k/k c . 



MFT do not depend on the frequency realizations) . (The 
perturbation theory Eq. I|16f> agreed with the frequency 
distribution approximation and was left out for clarity.) 
The error bars correspond to one standard deviation of 
the sample of 10 realizations. We note that the larger 
error bars occur after the transition. When the values of 
the order parameter are averaged over 10 realizations of 
the network and the frequencies, the results show very 
good agreement with the frequency distribution approx- 
imation and the directed mean field theory. 



In order to study how well our theory describes single 
realizations, we show in Fig.f^b) the order parameter r 2 
obtained from numerical solution of Eq. for a partic- 
ular realization of the network and frequencies (boxes), 
the time averaged theory (short dashed line) , and the fre- 
quency distribution approximation (solid line) as a func- 
tion of k/k c . As can be observed from the figure, in 
contrast with the time averaged theory, the frequency 
distribution approximation deviates from the numerical 
solution (boxes) by a small but noticeable amount. This 
behavior is observed for the other realizations as well. 
We note that the FDA and MFT results are virtually 
identical for all 10 realizations. On the other hand, the 
TAT and the results from direct numerical solution of 
Eq. show dependence on the realization. Since the 
FDA and MFT incorporate the realizations of the con- 
nections A nm , but not the frequencies, we interpret the 
observed realization dependence of the TAT and the di- 
rect solutions of Eq. Q as indicating that the latter de- 
pendence is due primarily to fluctuations in the realiza- 
tions of the frequencies rather than to fluctuations in the 
realizations of A„ m . 



Note that for our example N — 1500 and s = 2/15 im- 
plies that on average we have d m « d out « 200. Thus for 
comparison purposes, we generated an undirected net- 
work as follows: starting with N — 1500 nodes, we join 
pairs of nodes with undirected links in such a way that 
all nodes have = d° u * = 200. This is accomplished 
by using the configuration model described in Sec. IV of 
Ref. [l5|. The resulting network is described by a sym- 
metric adjacency matrix A. The results for this network 
are similar to those shown in the previous example. This 
suggests that the asymmetric network in the previous ex- 
ample can be considered (in a statistical sense) as sym- 
metric. 



and s — 2/15, Fig. fUa) shows the average of the order 
parameter r 2 obtained from numerical solution of Eq. l[TJl 
averaged over 10 realizations of the network and frequen- 
cies (triangles), the frequency distribution approximation 
(FDA, solid line), and the mean field theory (MFT, long 
dashed line) as a function of fc/fc c , where the results for 
the FDA and the MFT are averaged over the 10 net- 
work realizations (note, however, that the FDA and the 



In summary, for the random asymmetric network in 
example (i) and for the symmetric network described in 
the previous paragraph (not shown), all the approxima- 
tions work satisfactorily: single realizations are described 
by the time averaged theory, and the average over many 
realizations is described by the frequency distribution ap- 
proximation or the directed mean field theory. 
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Example (ii), A Strongly Asymmetric Network 
with A n n i > 



Now we consider a network in which the asymmetry 
has a more pronounced effect [example (ii)] . We consider 
directed networks defined in the following way. Using the 
configuration model as above, we first randomly generate 
an undirected network with N — 1500 nodes and 400 con- 
nections to each node, obtaining a symmetric adjacency 
matrix A' with entries or 1. We construct directed 
networks from this undirected network as follows. From 
the symmetric matrix A' , l's above the diagonal are in- 
dependently converted into 0's with probability 1 — p, 
generating by this process an asymmetric adjacency ma- 
trix A. (Imagining that the nodes are arranged in order 
of ascending n along a line, connections pointing in the 
direction of increasing n are randomly removed. This 
could model, for example, oscillators which are coupled 
chemically along the flow of some medium, or flashing 
fireflies that are looking mostly in one direction.) We 
will consider a rather low value of p, p = 0.1, in order to 
obtain a network with a strong asymmetry. 

In Fig we compare our approximations against the 
values of the order parameter obtained from numerical 
solution of Eq. Q as a function of k/k c for a network con- 
structed as described above where k c is given by Eq. (fT5|l . 
In Fig. Ola) we show the average of the order parameter 
r 2 [defined by Eq. (J3J] versus k/k c obtained from nu- 
merical solution of Eq. over 10 realizations of the 
network and frequencies (triangles), the frequency dis- 
tribution approximation (solid line), the directed mean 
field theory Eq. (|28l) (long dashed line) and the directed 
perturbation theory Eq. i|26|) (dotted-dashed line). The 
frequency distribution approximation captures, as in the 
undirected case, the values of the average of the order 
parameter obtained from numerical solution of Eq. IjTJl. 
The directed perturbation theory gives a good approxi- 
mation for small values of k close to k c , as expected. On 
the other hand, the directed mean field theory predicts 
a transition point which is smaller than the one actu- 
ally observed. We note that for this network solutions 
of Eq. yield substantial rms deviation of individual 
realizations [the error bars in Fig. OJa)] for all k > k c . 

Now we consider a single realization. In Fig. ^b) we 
show the order parameter r 2 obtained from numerical 
solution of Eq. for a particular realization of the net- 
work and frequencies (boxes), the time averaged theory 
(short dashed line) and the frequency distribution ap- 
proximation (solid line) as a function of k/k c . The time 
averaged theory tracks the value of the order parameter 
for this particular realization. This is also observed for 
the other realizations. 

As an indication of why the directed mean field theory 
gives a smaller transition point than that given by k c in 
Eq. (|15fl . we note that in the limiting case, p — > 0, all the 
elements above and in the diagonal of A are 0, so that 
A = and k c — > oo. However, the directed mean field 
theory predicts a transition at the finite value k m f = 
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FIG. 2: (a) Average of the order parameter r 2 obtained from 
numerical solution of Eq. Q over 10 realizations of the net- 
work and frequencies with p = 0.1 (triangles), from the fre- 
quency distribution approximation (solid line), from the di- 
rected mean field theory (long dashed line), and from the 
directed perturbation theory (dotted-dashed line) as a func- 
tion of k/k c . (b) Order parameter r 2 obtained from numerical 
solution of Eq. Q for a particular realization of the network 
and frequencies (boxes), from the time averaged theory (short 
dashed line) and from the frequency distribution approxima- 
tion (solid line) as a function of k/k c . 



ko(d in )/((d in d out )). 



C. Examples of Networks with Negative Coupling 

Now we consider examples in which there are negative 
connections, i.e., some of the entries of the adjacency 
matrix are negative, A nm < 0. In our next example, 
we construct first an undirected network with ./V = 1500 
nodes and 400 connections per node. We then set A nm = 
if n and m are not connected, and if they are we set 
A nm to 1 with probability q and to —1 with probability 
l-q. 

First we consider the case q = 2/3, so that one third 
of the connections are negative [example (iii)] . In Fig. [31 
we compare the numerical solution of Eq. Q with our 
theoretical approximations in Eqs. I|32|) and (|36|l for ten 
realizations of the network and frequencies. We show the 
average of the order parameter r 2 over 10 realizations of 
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FIG. 3: Average of the order parameter r 2 obtained from nu- 
merical solution of Eq. Q over 10 realizations of the network 
with q — 2/3 and frequencies (triangles with thin error bars), 
average of the time averaged theory (solid line with oval er- 
ror bars), and frequency distribution approximation (dashed 
line) as a function of k/k c . The horizontal line represents 
the value of the order parameter if the oscillators were phase 
locked {On = 6m for all m and n). 



the network (triangles with thin error bars), the average 
of the TAT [Eq. I(32|) , solid line with oval error bars] , and 
the average of the FDA [Eq. ((36(1 , dashed line] . The er- 
ror bar widths represent one standard deviation of the 
sample of 10 realizations. As in the previous examples, 
the FDA did not show noticeable variations for different 
realizations of the network. We observe that the order 
parameter computed from our theory yields a slightly 
larger value than that obtained from the numerical solu- 
tion of Eq. , but in general both the transition point 
and the behavior of the order parameter are described 
satisfactorily by the theory. 

In this case, the phases ip n obtained from numerical 
solution of Eq. (|32[) do not depend on n, i.e., ip n = ip m for 
all n, to. This can be understood on the basis that there 
are not enough negative coupling terms to make the right 
hand side of Eq. 1(12(1 negative, so that a solution exists 
in which all the phases ip n are equal. As mentioned in 
Sec. II VI the difference in the phases in Eq. ((32(1 prevents 
the right hand side of Eq. 112|) from becoming negative in 
the presence of negative connections. As a confirmation 
of this we note that as k — > oo the order parameter r 
appears to approach 1/3 (the dotted-dashed horizontal 
line in Fig. which corresponds to (tp n — ip m ) — > in 
Eq. 1(35(1 for q = 2/3. The fact that both the phases ipn 
and 8 n do not depend on n as k — > oo is consistent with 
Eq. J7J. 

In order to consider a case in which the effect of the 
negative connections is more extreme, we consider a net- 
work constructed as described above with with q = 0.54 
[example (iv)]. In Fig. 0] we compare the numerical so- 
lution of Eq. Q with our theoretical approximations in 
Eqs. (|32|) and ((36(1 for ten realizations of the network 
and frequencies. We show the average of the order pa- 



rameter r 2 over 10 realizations of the network (triangles 
with thin error bars), the average of the FDA [Eq. 136( 1. 
dashed line with thin error bars] and the average of the 
TAT [Eq. J22J), solid line with oval error bars]. When 
numerically solving Eq. (|32|l by iteration of Eq. 1(33(1 . on 
some occasions a period two orbit was found instead of 
the desired fixed point. If we denote the left hand side 
of Eq. (|33|l by and the right hand side by f{z^), we 
found that convergence to a fixed point was facilitated 
by replacing the right hand side by [z 3 n + f{z 3 n )]/2 and 
finding the fixed points of this modified system. 
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FIG. 4: Average of the order parameter r 2 obtained from nu- 
merical solution of Eq. £Q| over 10 realizations of the network 
with q — 0.54 and frequencies (triangles) and average of the 
TAT (solid line) as a function of k/k c . Note the different scale 
in the horizontal axis as compared with the previous figures. 
The horizontal dotted-dashed line represents the value of the 
order parameter if the oscillators were phase locked (6 n = 8 m 
for all m and n). 

In this example, at low coupling strengths [roughly 
k/k c < 4, where k c is computed from Eq. (|37|) ] the order 
parameter computed from numerical solution of Eq. Q is 
smaller than that obtained from the TAT and FDA. As k 
increases, however, the TAT and FDA theories captures 
the asymptotic value of the order parameter r. We note 
that in this case the asymptotic value is larger than that 
corresponding to phase locking [i.e., the one obtained by 
setting tfj n = in Eq. ®, r w 0.54-0.46 = 0.08], which 
we indicate by a horizontal dotted-dashed line in Fig^] 
and much smaller than r — 1, the value corresponding to 
no frustration [i.e., tj) n — ip m — for A nm > and ir for 
A nm < in Eq. 1)35(1] . The small scale of the horizontal 
axis is due to the fact that we are plotting r 2 , and to our 
definition of the order parameter which assigns a value 
of 1 to a non frustrated configuration. The small value 
of the order parameter indicates a strong frustration. 

We note that in this example, in contrast with the ex- 
amples discussed so far, there is variation in the values of 
the order parameter predicted by the FDA for different 
realizations of the network. This indicates that, as the 
expected value of the coupling strengths A nm becomes 
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small (i.e., \q — 1/2 1 small), fluctuations due to the re- 
alization of the network become noticeable. Although 
the values predicted by the FDA and TAT depend on 
the realization of the network and frequencies, we note 
for k/k c > 6 that these values track the values observed 
for the numerical simulations of the corresponding real- 
ization. As an illustration of this, we plot in Fig. |SJa) 
the values of r 2 obtained from the TAT (triangles) and 
in Fig. Efb) the values of r 2 obtained from the FDA 
(boxes) versus the value obtained from numerical solu- 
tion of Eq. (JIJ for k/k c = 8. Besides a small positive bias 
in the FDA, the theories track the spread in the results 
of the numerical solution for different realizations. Some 
bias in the FDA is not surprising, because we averaged 
the right hand side of the nonlinear equation (|12|l for the 
TAT in order to get Eq. (T^J for the FDA. Nonetheless, 
the bias is extremely small in most of our examples. 
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FIG. 5: Order parameter r 2 obtained from (a) the TAT (tri- 
angles) and from (b) the FDA (boxes) versus the value ob- 
tained from numerical solution of Eq. Q for k/k c = 8. The 
solid line is the identity. Besides a small positive bias in the 
FDA, the theories track the spread in the results of the nu- 
merical solution for different realizations. 



The behavior observed in Fig. 0] at k/k c < 4 can be 
interpreted as a shift in the transition point to a larger 



value of the coupling strength, and is reminiscent of what 
occurs when the time fluctuations [kh n (i] in Eq. (JSJ] ne- 
glected in Eq. © have an appreciable effect. We be- 
lieve that the time fluctuations have a more pronounced 
effect as the number of negative connections becomes 
comparable to the number of positive connections (i.e., 
as | q — 1/2 1 becomes small) because the critical coupling 
strength k c becomes large (roughly k c ^ \q— 1 /2 1 1 ) . In 
particular, with positive connections, the condition for 
neglecting kh n (t) was that the number of connections 
to each node was large. In contrast, for the present case, 
the analogous statement would be that \q— 1/2 1 times the 
number of connections is large, which is much less well 
satisfied, \q - 1/2|400 = 0.04 x 400 = 16. The extreme 
case of zero mean coupling has already been studied nu- 
merically by Daidoi^, who found that in this case the os- 
cillators lock in the sense that their average frequenc y is 
the same, but their phases diffuse. As argued in Ref. [T3, 
such fluctuations have the effect of shifting the transi- 
tion to larger values of the coupling strength. It would 
be interesting to carry on simulations in networks with 
a much larger number of connections per node, as the 
effect of fluctuations would likely be reduced. 

We also considered a case in which the adjacency ma- 
trix is asymmetric and has mixed positive/negative con- 
nections. For N = 1500 nodes, we constructed an adja- 
cency matrix by setting its nondiagonal entries to 1 , — 1 , 
and with probability 8/45, 4,45, and 11/15, respec- 
tively. The latter probability yields an expected number 
of connections of 400. Our theories work satisfactorily 
in this case, and, since the results are similar to those 
in Fig. |31 we do not show them. In this case there is no 
guarantee that there is a real eigenvalue [as needed for 
estimating the critical coupling strength in Eq. I|15(l ]. or 
that the largest real eigenvalue (if there is one) has the 
largest real part. Numerically, we find that for matrices 
constructed as in this example there is a real positive 
eigenvalue and that, furthermore, it is well separated 
from the largest real part of the remaining eigenvalues 
(see Fig. EJ). We also find this for other values of q pro- 
vided I q — 1/2 1 is not too small. We provide a discussion 
of this issue and show the spectrum of the adjacency ma- 
trix in Appendix lAl 



VI. DISCUSSION 

In this paper, we have considered interacting phase 
oscillators [Eq. QJ] connected by directed networks and 
networks with mixed positive /negative connections. 

The previous theory of Ref. [jjl given by Eq. Q12[l (the 
time averaged theory, TAT) can still be applied for asym- 
metric networks with purely positive coupling and was 
found to give good predictions, applicable to individual 
asymmetric random realizations [Figs. ^b),[2Ib)]. The 
previous theory given by Eq. I|13|) (the frequency distri- 
bution approximation, FDA) can also still be applied 
for asymmetric networks with purely positive coupling 
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and was found to give good predictions applicable to the 
ensemble average behavior of asymmetric network real- 
izations [Figs, da), EJa)]. The perturbative theory for 
the FDA was generalized to account for directed net- 
works [Eqs. J2j and as was the previous undi- 
rected network mean field theory, MFT (generalized from 

Eqs. (JTHJ-lJiBj to Ec l s - In our example (ii), 

which had a very strong asymmetry, we found that our 
directed FDA perturbation theory [Eqs. and (j?7jl ] 
gave a good description of synchronization, but that the 
directed mean field approximation gave a transition to 
synchronization at a coupling substantially below that 
observed. In contrast, for example (i), in which the cou- 
pling matrices were individually asymmetric but their 
ensemble average was symmetric, the mean field theory 
(and all the other theories in Sec. 1111(1 gave good results. 

For the case of mixed positive/negative couplings 
we presented a generalization of the TAT and FDA, 
Eqs. (J22J)-lj2ZI). We tested these results on two exam- 
ples, example (iii) in which a fraction 1 — q = 1/3 of 
the couplings were negative, and example (iv) in which a 
fraction 1 — q = 0.46 of the couplings were negative. For 
example (iii) we found that iteration of Eq. (|33fl converges 
to a fixed point with ip n — ip m = 0, and thus the result 
is similar to the case where all connections are positive. 
In example (iv), the result of iteration of Eq. H33[) yields 
nontrivial values for the phases ip n . In this case we found 
good agreement between the solutions of Q and the the- 
ory for the order parameter for k/k c large (k/k c > 4), but 
that for smaller k/k c (k/k c < 4), although yielding qual- 
itatively similar behavior to that observed (Fig. , the 
theory overestimates the order parameter. Analogous to 
similar observations for symmetric networks with only 
positive coupling^, we speculate (Sec. 0) that this is 
a finite size effect associated with the fact that the ef- 
fective number of connections given in this example by 
\q— 1/2 1 400 = 16 is not sufficiently large to justify neglect 
of kh n (t) in Eq. © 

In order to isolate the effect of the asymmetry and the 
negative connections, we considered networks in which 
the degree distribution is very narrow. The combined 
effect of these factors with different heterogeneous de- 
gree distributions (e.g., scale free networks^) and with 
correlations in the network (in particular, degree-degree 
correlations) is still open to investigation. 

In practice, one could be interested in networks in 
which the asymmetry in the connections is strongly cor- 
related with the sign of the coupling (in analogy to some 
models in neuroscience^) . Although we did not study 
such a case here, we believe our theory provides a good 
starting point to study the emergence of synchronization 
in these kind of structured complex networks. 



APPENDIX A 

In this Appendix we discuss the characteristics of the 
spectrum of the adjacency matrices considered in our ex- 



amples. Although we will focus here on asymmetric ma- 
trices, a similar argument works for symmetric matrices. 
The matrices we consider are relatively sparse, with the 
position of the nonzero entries being chosen randomly 
(e.g., in the symmetric case, the position of the nonzero 
entries is chosen when constructing the network using 
the configuration model) , and their values being also de- 
termined randomly from a given probability distribution 
(e.g., 1 with probability q and —1 with probability 1 — q). 
Our interest is focused on the gap between the largest 
real eigenvalue (if there is one) and the largest real part 
of the other eigenvalues. In Ref. the spectrum of cer- 
tain large sparse matrices with average eigenvalue and 
row sum Yl m =i ^nm = 1 was described and a heuristic 
analytical approach was proposed. Using results for ma- 
trices with zero mean Gaussian random entries^, Ref. 
I23I predicts that the spectrum of the non-Gaussian ran- 
dom matrices they consider consists of a trivial eigen- 
value A = 1 with the remaining eigenvalues distributed 
uniformly in a circle centered at the origin of the complex 
plane with radius 

p = VNa, (Al) 
where a 1 is the variance of the entries of the matrix. We 
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FIG. 6: Complex eigenvalues A (dots) of a 1500 x 1500 ran- 
dom matrix whose off diagonal entries are 1, —1 or with 
probabilities 8/45, 4/45, and 11/15, respectively. One eigen- 
value is located at A = 136.2, while the other 1499 eigenvalues 
uniformly fill a circle of radius p centered at the origin of the 
complex A plane. Note that p « 19.8 is substantially less than 
136.2. Comparing with the theory in the Appendix, Eq. (IA2II 
yields a prediction of 133.3 for the maximum real eigenvalue 
while Eq. 1A11 predicts 19.7 for p. These are in excellent 
agreement with our numerically determined values. 



find that this approach also succeeds in describing the 
spectrum of the matrices in our examples. In our case, 
the diagonal entries are 0, so that the average eigenvalue 
is also as in Ref. |23- We find that there is always a 
largest real eigenvalue approximately given by the mean 
field value 

\ = {d 2 )/{d) (A2) 

(see Refs. [T2j2j), where d n = ^ =1 i„ ro and (d 2 ) = 
J2n=i which in the case considered in Ref. 2^ reduces 
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to A = 1 . We also numerically confirm that the remaining 
eigenvalues arc uniformly distributed in a circle of radius 
p as described in Ref. |23 This is illustrated in Fig. H3 

Thus for N 3> 1 if A > p there is a gap of size A — 
p between the largest real eigenvalue and real part of 
the rest of the eigenvalue spectrum. Using Eqs. (|A1|> 
and \A2\ it can be shown that, for networks with large 
enough number of connections per node or with enough 



positive (or negative) bias in the coupling strength, there 
is a wide separation between the largest eigenvalue and 
the largest real part of the remaining eigenvectors. For 
symmetric matrices, similar results apply (i.e, the bulk 
of the spectrum of the matrix A can be approximately 
obtained as described above using Wigner's semicircle 
law). 
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